Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix roundoff issue in SUNDIALS evolve() #4148

Merged
merged 5 commits into from
Sep 17, 2024

Conversation

ajnonaka
Copy link
Contributor

Summary

Fix a roundoff issue in the SUNDIALS evolve() function that caused the integrator to run a second time with a very tiny dt in a given time step.

Additional background

Checklist

The proposed changes:

  • [X ] fix a bug or incorrect behavior in AMReX
  • add new capabilities to AMReX
  • changes answers in the test suite to more than roundoff level
  • are likely to significantly affect the results of downstream AMReX users
  • include documentation in the code and/or rst files, if appropriate

@ajnonaka ajnonaka changed the title Fix foundoff issues in SUNDIALS evolve() Fix roundoff issues in SUNDIALS evolve() Sep 12, 2024
@ajnonaka ajnonaka changed the title Fix roundoff issues in SUNDIALS evolve() Fix roundoff issue in SUNDIALS evolve() Sep 12, 2024
@@ -251,7 +251,8 @@ public:
for (int step_number = 0; step_number < BaseT::max_steps && !stop; ++step_number)
{
// Adjust step size to reach output time
if (time_out - time_current < dt) {
// protect against roundoff
if (std::abs(time_out - time_current - dt) <= 1.e5 * std::numeric_limits<amrex::Real>::epsilon() * std::abs(dt)) {
Copy link
Member

@WeiqunZhang WeiqunZhang Sep 12, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What if Real is float?

We could also do amrex::almostEqual(time_current+dt, time_out). almostEqual has an optional argument that can be used to fine tune what is considered almost equal.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That sounds good, I'll test it.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@WeiqunZhang how is this?

@WeiqunZhang
Copy link
Member

Any reasons why std::abs is needed? Also the last argument should be an int (e.g., 10, not double like 10..

@ajnonaka
Copy link
Contributor Author

@WeiqunZhang fixed

@WeiqunZhang WeiqunZhang merged commit 73716d9 into AMReX-Codes:development Sep 17, 2024
59 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants